W pliku banknotes.csv znadjują się dane opisujące obrazy
banknotów. Dane powstały poprzez transformatę falową (wavelett
transform) zastosowaną do obrazów w skali szarości rozmiaru 400x400
pikseli. Po zastosowaniu transformaty wyliczono cztery charakterystyki
liczbowe obrazu - wariancję, skośność, kurtozę oraz entropię.
Za pomocą modelu regresji logistycznej sprawdź czy za pomocą tej metody jesteśmy w stanie dobrze odróżnić banknoty prawdziwe od fałszywych.
library(tidyverse)
## Warning: pakiet 'tidyverse' został zbudowany w wersji R 4.2.3
## Warning: pakiet 'ggplot2' został zbudowany w wersji R 4.2.3
## Warning: pakiet 'readr' został zbudowany w wersji R 4.2.3
## Warning: pakiet 'forcats' został zbudowany w wersji R 4.2.3
## Warning: pakiet 'lubridate' został zbudowany w wersji R 4.2.3
library(caret)
## Warning: pakiet 'caret' został zbudowany w wersji R 4.2.3
library(ggfortify)
## Warning: pakiet 'ggfortify' został zbudowany w wersji R 4.2.3
library(knitr)
library(kableExtra)
## Warning: pakiet 'kableExtra' został zbudowany w wersji R 4.2.3
library(plotly)
## Warning: pakiet 'plotly' został zbudowany w wersji R 4.2.3
bank <- readr::read_csv('C:/Users/PC/Downloads/banknote.csv', col_names = FALSE)
Dane zapisane są w następującej postaci:
X1 - wariancja,X2 - skośność,X3 - kurtoza,X4 - entropia,X5 - 1 prawdziwe, 0 fałszywe.bank %>% count(X5)
Widzimy, że zbiór danych banknote nie jest zbilansowany.
Posiada więcej fałszywych banknotów niż prawdziwych.
Przy użyciu histogramów oraz testu Shapiro-Wilka sprawdzamy rozkłady naszych zmiennych.
wykres <- plot_ly(alpha = 0.7)
wykres <- wykres %>% add_histogram(x = bank$X1, name = "X1")
wykres <- wykres %>% add_histogram(x = bank$X2, name = "X2")
wykres <- wykres %>% add_histogram(x = bank$X3, name = "X3")
wykres <- wykres %>% add_histogram(x = bank$X4, name = "X4")
wykres <- wykres %>% layout(barmode = "stack")
wykres
shapiro.test(bank$X1)
##
## Shapiro-Wilk normality test
##
## data: bank$X1
## W = 0.982, p-value = 4.686e-12
shapiro.test(bank$X2)
##
## Shapiro-Wilk normality test
##
## data: bank$X2
## W = 0.97463, p-value = 8.224e-15
shapiro.test(bank$X3)
##
## Shapiro-Wilk normality test
##
## data: bank$X3
## W = 0.92709, p-value < 2.2e-16
shapiro.test(bank$X4)
##
## Shapiro-Wilk normality test
##
## data: bank$X4
## W = 0.91488, p-value < 2.2e-16
Z testu Shapiro-Wilka \(p-value < 0,05\), zatem odrzucamy \(H_0\) i stwierdzamy, że wariancja, skośność, kurtoza oraz entropia nie mają rozkładu normalnego.
regresja_logistyczna <- glm(X5~X1+X2+X3+X4,data=bank,family = "binomial")
## Warning: glm.fit: dopasowane prawdopodobieństwa numerycznie okazały się być 0
## lub 1
summary(regresja_logistyczna)
##
## Call:
## glm(formula = X5 ~ X1 + X2 + X3 + X4, family = "binomial", data = bank)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.70001 0.00000 0.00000 0.00029 2.24614
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 7.3218 1.5589 4.697 2.64e-06 ***
## X1 -7.8593 1.7383 -4.521 6.15e-06 ***
## X2 -4.1910 0.9041 -4.635 3.56e-06 ***
## X3 -5.2874 1.1612 -4.553 5.28e-06 ***
## X4 -0.6053 0.3307 -1.830 0.0672 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 1885.122 on 1371 degrees of freedom
## Residual deviance: 49.891 on 1367 degrees of freedom
## AIC: 59.891
##
## Number of Fisher Scoring iterations: 12
Estimate: szacunki przechwycenia \((\beta_0)\) i współczynnika beta związanego
z każdą zmienną przewidującąStd.Error: błąd standardowy oszacowań współczynników.
Przedstawia on dokładność współczynników. Im większy błąd standardowy,
tym mniejsza pewność co do oszacowania.wartość z: statystyka Z, która jest estymacją
współczynnika (kolumna 2) podzieloną przez błąd standardowy estymacji
(kolumna 3)Pr(>|z|): Wartość p odpowiadająca statystyce z. Im
mniejsza wartość p, tym bardziej znaczący jest szacunek.Różnica pomiędzy odchyleniem zerowym a odchyleniem resztowym mówi
nam, że model jest dobrze dopasowany. Większa różnica jest lepsza dla
modelu. Null deviance to wartość, gdy w równaniu mamy tylko
przechwyt bez żadnych zmiennych, a Residual deviance to
wartość, gdy bierzemy pod uwagę wszystkie zmienne. Ma to sens, aby uznać
model za dobry, jeśli ta różnica jest wystarczająco duża.
W naszym przypadku wartość funkcji logitowej \(ln\bigg(\frac{p}{1-p}\bigg)\) z autentyczności banknotów dla wariancji, skośności, kurtozy oraz entropii, wynosi \(7,3218\).
dec_bond <- regresja_logistyczna$coefficients
logistic2_slope = -dec_bond[2]/dec_bond[3]
logistic2_intercept = -dec_bond[1]/dec_bond[3]
ggplot(bank, aes(x = X1, y = X2, color = as.factor(X5))) +
geom_point() +
geom_abline(slope = logistic2_slope, intercept = logistic2_intercept, color = 'blue', linetype = 'dashed') +
labs(color = 'X5')
ggplot(bank, aes(x = X1, y = X3, color = as.factor(X5))) +
geom_point() +
geom_abline(slope = logistic2_slope, intercept = logistic2_intercept, color = 'blue', linetype = 'dashed') +
labs(color = 'X5')
ggplot(bank, aes(x = X1, y = X4, color = as.factor(X5))) +
geom_point() +
geom_abline(slope = logistic2_slope, intercept = logistic2_intercept, color = 'blue', linetype = 'dashed') +
labs(color = 'X5')
Wykresy przedstawiają hiperpłaszyznę decyzyjną naszego modelu.
logit_prediction <- predict(regresja_logistyczna, bank, type = 'response')
head(logit_prediction)
## 1 2 3 4 5 6
## 2.220446e-16 2.220446e-16 2.185822e-10 2.220446e-16 4.579103e-01 2.220446e-16
logistic2_predictions <- ifelse(logit_prediction > 0.5, 1, 0)
head(logistic2_predictions)
## 1 2 3 4 5 6
## 0 0 0 0 0 0
library(caret)
caret::confusionMatrix(data = as.factor(logistic2_predictions), reference = as.factor(bank$X5))
## Confusion Matrix and Statistics
##
## Reference
## Prediction 0 1
## 0 757 6
## 1 5 604
##
## Accuracy : 0.992
## 95% CI : (0.9857, 0.996)
## No Information Rate : 0.5554
## P-Value [Acc > NIR] : <2e-16
##
## Kappa : 0.9838
##
## Mcnemar's Test P-Value : 1
##
## Sensitivity : 0.9934
## Specificity : 0.9902
## Pos Pred Value : 0.9921
## Neg Pred Value : 0.9918
## Prevalence : 0.5554
## Detection Rate : 0.5517
## Detection Prevalence : 0.5561
## Balanced Accuracy : 0.9918
##
## 'Positive' Class : 0
##
table <- data.frame(confusionMatrix(as.factor(logistic2_predictions), as.factor(bank$X5))$table)
plotTable <- table %>%
mutate(goodbad = ifelse(table$Prediction == table$Reference, "good", "bad")) %>%
group_by(Reference) %>%
mutate(prop = Freq/sum(Freq))
ggplot(data = plotTable, mapping = aes(x = Reference, y = Prediction, fill = goodbad, alpha = prop)) +
geom_tile() +
geom_text(aes(label = Freq), vjust = .5, fontface = "bold", alpha = 1) +
scale_fill_manual(values = c(good = "green", bad = "red")) +
theme_bw() +
labs(title = "Macierz modelu regresja_logistyczna") +
theme(plot.title = element_text(hjust = 0.5)) +
xlim(rev(levels(table$Reference)))
Przy pomocy biblioteki caret stworzyliśmy macierz
pomyłek. Jej charakterystyki liczbowe przedstawiają się w następujący
sposób:
Dokładność to liczba próbek poprawnie sklasyfikowanych spośród wszystkich próbek obecnych w zbiorze testowym:
\(Accuracy = \frac{TP+TN}{TP+FP+TN+FN}\), u nas wynosi \(0,992\).
\(NIR > 0,5\), zatem nasz model nie jest gorszy od rzutu monetą.
Współczynnik Kappa może być użyty do oceny dokładności klasyfikacji. Przyjmuje wartości \([-1,1]\) wartość większa niż 0 wskazuje, że klasyfikacja jest znacznie lepsza niż losowa.
Czułość to liczba próbek rzeczywiście należących do klasy pozytywnej spośród wszystkich próbek, które zostały przewidziane przez model jako należące do klasy pozytywnej.
\(Sensitivity = \frac{TP}{TP+FN}\).
Swoistość to liczba próbek przewidywanych prawidłowo jako należące do klasy negatywnej spośród wszystkich próbek w zbiorze danych, które rzeczywiście należą do klasy negatywnej.
\(Sensitivity = \frac{TN}{TN+FP}\).
library(caret)
train_test_split <- createDataPartition(bank$X5, list = FALSE, p=0.75)
bank_train <- bank[train_test_split,]
bank_test <- bank[-train_test_split,]
cat(dim(bank_train),dim(bank_test))
## 1029 5 343 5
table(bank$X5)
##
## 0 1
## 762 610
Wszystkie banknoty prawdziwe i fałszywe dostępne.
table(bank_train$X5)
##
## 0 1
## 568 461
Banknoty treningowe, tylko \(\frac{3}{4}\) banknotów. (p=0.75)
table(bank_test$X5)
##
## 0 1
## 194 149
Te banknoty, które odjeliśmy od wszystkich banknotów by otrzymać banknoty treningowe.
kable(summary(bank_train), escape = F) %>%
kable_styling(bootstrap_options = c("striped", "hover", "condensed"))
| X1 | X2 | X3 | X4 | X5 | |
|---|---|---|---|---|---|
| Min. :-7.0421 | Min. :-13.773 | Min. :-5.2861 | Min. :-8.5482 | Min. :0.000 | |
| 1st Qu.:-1.8439 | 1st Qu.: -2.125 | 1st Qu.:-1.4801 | 1st Qu.:-2.1888 | 1st Qu.:0.000 | |
| Median : 0.3798 | Median : 2.084 | Median : 0.6766 | Median :-0.5724 | Median :0.000 | |
| Mean : 0.4016 | Mean : 1.777 | Mean : 1.4926 | Mean :-1.1548 | Mean :0.448 | |
| 3rd Qu.: 2.8237 | 3rd Qu.: 6.777 | 3rd Qu.: 3.4356 | 3rd Qu.: 0.3948 | 3rd Qu.:1.000 | |
| Max. : 6.5633 | Max. : 12.952 | Max. :17.9274 | Max. : 2.4495 | Max. :1.000 |
bank_train_model <- glm(X5 ~ X1+X2+X3+X4 , data = bank_train, family='binomial')
## Warning: glm.fit: dopasowane prawdopodobieństwa numerycznie okazały się być 0
## lub 1
summary(bank_train_model)
##
## Call:
## glm(formula = X5 ~ X1 + X2 + X3 + X4, family = "binomial", data = bank_train)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.84582 0.00000 0.00000 0.00047 2.04983
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 6.8728 1.7249 3.984 6.76e-05 ***
## X1 -7.2392 1.9281 -3.755 0.000174 ***
## X2 -4.0288 1.0298 -3.912 9.15e-05 ***
## X3 -4.9895 1.3060 -3.820 0.000133 ***
## X4 -0.6681 0.3912 -1.708 0.087669 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 1415.35 on 1028 degrees of freedom
## Residual deviance: 35.17 on 1024 degrees of freedom
## AIC: 45.17
##
## Number of Fisher Scoring iterations: 12
logit_prediction_bank_train_model <- predict(bank_train_model, bank_train, type = 'response')
head(logit_prediction_bank_train_model)
## 1 2 3 4 5 6
## 2.220446e-16 2.220446e-16 1.762616e-09 2.220446e-16 5.715952e-01 2.220446e-16
logistic2_predictions_bank_train_model <- ifelse(logit_prediction_bank_train_model > 0.5, 1, 0)
head(logistic2_predictions_bank_train_model)
## 1 2 3 4 5 6
## 0 0 0 0 1 0
library(caret)
caret::confusionMatrix(data = as.factor(logistic2_predictions_bank_train_model), reference = as.factor(bank_train$X5))
## Confusion Matrix and Statistics
##
## Reference
## Prediction 0 1
## 0 562 5
## 1 6 456
##
## Accuracy : 0.9893
## 95% CI : (0.981, 0.9947)
## No Information Rate : 0.552
## P-Value [Acc > NIR] : <2e-16
##
## Kappa : 0.9784
##
## Mcnemar's Test P-Value : 1
##
## Sensitivity : 0.9894
## Specificity : 0.9892
## Pos Pred Value : 0.9912
## Neg Pred Value : 0.9870
## Prevalence : 0.5520
## Detection Rate : 0.5462
## Detection Prevalence : 0.5510
## Balanced Accuracy : 0.9893
##
## 'Positive' Class : 0
##
table <- data.frame(confusionMatrix(as.factor(logistic2_predictions_bank_train_model), as.factor(bank_train$X5))$table)
plotTable <- table %>%
mutate(goodbad = ifelse(table$Prediction == table$Reference, "good", "bad")) %>%
group_by(Reference) %>%
mutate(prop = Freq/sum(Freq))
ggplot(data = plotTable, mapping = aes(x = Reference, y = Prediction, fill = goodbad, alpha = prop)) +
geom_tile() +
geom_text(aes(label = Freq), vjust = .5, fontface = "bold", alpha = 1) +
scale_fill_manual(values = c(good = "green", bad = "red")) +
theme_bw() +
labs(title = "Macierz modelu bank_train") +
theme(plot.title = element_text(hjust = 0.5)) +
xlim(rev(levels(table$Reference)))
fmpreds <- predict(bank_train_model, bank_train, type = 'response')
fmpreds_classes <- ifelse(fmpreds > 0.5, 1, 0)
baseline_cm <- confusionMatrix(factor(fmpreds_classes), factor(bank_train$X5))
baseline_cm
## Confusion Matrix and Statistics
##
## Reference
## Prediction 0 1
## 0 562 5
## 1 6 456
##
## Accuracy : 0.9893
## 95% CI : (0.981, 0.9947)
## No Information Rate : 0.552
## P-Value [Acc > NIR] : <2e-16
##
## Kappa : 0.9784
##
## Mcnemar's Test P-Value : 1
##
## Sensitivity : 0.9894
## Specificity : 0.9892
## Pos Pred Value : 0.9912
## Neg Pred Value : 0.9870
## Prevalence : 0.5520
## Detection Rate : 0.5462
## Detection Prevalence : 0.5510
## Balanced Accuracy : 0.9893
##
## 'Positive' Class : 0
##
Analiza krzywej charakterystyki operacyjnej odbiornika (ROC — Receiver operating characteristic) jest użytecznym narzędziem oceny dokładności predykcji modelu poprzez wykreślenie czułości wobec (1-swoistości) testu klasyfikacyjnego (przy progu zmieniającym się w całym zakresie wyników testu diagnostycznego).
roc_rl <- pROC::roc(response = bank_train$X5, predictor = fmpreds)
pROC::ggroc(roc_rl, legacy.axes = TRUE) + geom_abline(slope = 1, intercept = 0) +
labs(title = "Krzywa ROC bank_train")+
theme(plot.title = element_text(hjust = 0.5))
roc_rl$auc
## Area under the curve: 0.9999
Krzywa kształtem bardzo przypomina trójkąc co oznacza, że nasz model jest prawie “idealny”, a pole pod wykresem wynosi \(0,9998\).
Dzięki tej metodzie jesteśmy w stanie odróżnić prawdziwe banknoty od
fałszywych, ale zdarza się, że nasz model popełni błąd. Tworząc model
treningowy z \(\frac{3}{4}\) wszystkich
banknotów uzyskaliśmy lepsze współczynniki
No Information Rate, Sensitivity oraz
Specificity. Co za tym idzie nasz model jest bliższy
modelowi o idealnych parametrach.